Diffusion Tomography system and method using direct reconstruction of scattered radiation

ABSTRACT

A method for the direct reconstruction of an object from measurements of the transmitted intensity of diffusively scattered radiation effected by irradiating the object with a source of continuous wave radiation. The transmitted intensity is related to the diffusion coefficient by an integral operator. The image is directly reconstructed by executing a prescribed mathematical algorithm, as determined with reference to an integral operator, on the transmitted intensity of the diffusively scattered radiation.

FIELD OF THE INVENTION

This invention relates generally to a system, and concomitantmethodology, for generating an image of an object and, moreparticularly, to such system and methodology for which the image isdirectly reconstructed from measurements of scattered radiation detectedby irradiating the object with a continuous wave source.

BACKGROUND OF THE INVENTION

The inventive subject matter addresses the physical principles and theassociated mathematical formulations underlying the directreconstruction method for optical imaging in the multiple scatteringregime. The result is a methodology for the direct solution to the imagereconstruction problem.

Moreover, the method is generally applicable to imaging with any scalarwave in the diffusive multiple scattering regime and is not limited tooptical imaging. However, for the sake of elucidating the significantramifications of the present invention, it is most instructive to selectone area of application of the method so as to insure a measure ofdefiniteness and concreteness to the description. Accordingly, sincemany biological systems meet the physical requirements for theapplication of the principles of the present invention, especiallydiffusion tomography principles, the fundamental aspects of the presentinventive subject matter will be conveyed using medical imaging as anillustrative application of the method.

There have been three major developments in medical imaging over thepast 20 years that have aided in the diagnosis and treatment of numerousmedical conditions, particularly as applied to the human anatomy; thesedevelopments are: (1) the Computer-Assisted Tomography (CAT) scan; (2)the Magnetic Resonance Imaging (MRI); and (3) the Positron EmissionTomography (PET) scan.

With a CAT scanner, X-rays are transmitted through, for example, a humanbrain, and a computer uses X-rays detected external to the human head tocreate and display a series of images--basically cross-sections of thehuman brain. What is being imaged is the X-ray absorption function forunscattered, hard X-rays within the brain. CAT scans can detect, forinstance, strokes, tumors, and cancers. With an MRI device, a computerprocesses data from radio signals impinging on the brain to assemblelife-like, three-dimensional images. As with a CAT scan, suchmalformations as tumors, blood clots, and atrophied regions can bedetected. With a PET scanner, the positions of an injected radioactivesubstance are detected and imaged as the brain uses the substance. Whatis being imaged is the GAMMA ray source position. Each of these medicalimaging techniques has proved invaluable to the detection and diagnosingof many abnormal medical conditions. However, in many respects, none ofthe techniques is completely satisfactory for the reasons indicated inthe following discussion.

In establishing optimal design parameters for a medical imagingtechnique, the following four specifications are most important. Thespecifications are briefly presented in overview fashion before a moredetailed discussion is provided; moreover, the shortcomings of each ofthe conventional techniques are also outlined. First, it would bepreferable to use a non-ionizing source of radiation. Second, it wouldbe advantageous to achieve spatial resolution on the order of amillimeter to facilitate diagnosis. Third, it would be desirable toobtain metabolic information. And, fourth, it would be beneficial toproduce imaging information in essentially real-time (on the order ofone millisecond) so that moving picture-like images could be viewed.None of the three conventional imaging techniques is capable ofachieving all four specifications at once. For instance, a CAT scanneris capable of high resolution, but it uses ionizing radiation, it is notcapable of metabolic imaging, and its spatial resolution is borderlineacceptable. Also, while MRI does use non-ionizing radiation and hasacceptable resolution, MRI does not provide metabolic information and isnot particularly fast. Finally, a PET scanner does provide metabolicinformation, but PET uses ionizing radiation, is slow, and spatialresolution is also borderline acceptable. Moreover, the PET technique isinvasive due to the injected substance.

The four specifications are now considered in more detail. With respectto ionizing radiation, a good deal of controversy as to its effects onthe human body presently exists in the medical community. To ensure thatthe radiation levels are within what are now believed to be acceptablelimits, PET scans cannot be performed at close time intervals(oftentimes, it is necessary to wait at least 6 months between scans),and the dosage must be regulated. Moreover, PET is still a research toolbecause a cyclotron is needed to make the positron-emitting isotopes.Regarding spatial resolution, it is somewhat self-evident that diagnosiswill be difficult without the necessary granularity to differentiatedifferent structures as well as undesired conditions such as blood clotsor tumors. With regard to metabolic information, it would be desirable,for example, to make a spatial map of oxygen concentration in the humanhead, or a spatial map of glucose concentration in the brain. Theability to generate such maps can teach medical personnel about diseaseas well as normal functions. Unfortunately, CAT and MRI report densitymeasurements--electrons in an X-ray scanner or protons in MRI--and thereis not a great deal of contrast to ascertain metabolic information, thatis, it is virtually impossible to distinguish one chemical (such asglucose) from another. PET scanners have the ability to obtain metabolicinformation, which suggests the reason for the recent popularity of thistechnique. Finally, imaging is accomplished only after a substantialprocessing time, so real-time imaging is virtually impossible with theconventional techniques.

Because of the aforementioned difficulties and limitations, there hasbeen much current interest in the development of a technique forgenerating images of the distribution of diffusion coefficients ofliving tissue that satisfy the foregoing four desiderata. Accordingly, atechnique using low intensity photons would be safe. The techniqueshould be fast in that optical events occur within the range of 10nanoseconds--with this speed, numerous measurements could be completedand averaged to reduce measurement noise while still achieving the onemillisecond speed for real-time imaging. In addition, source anddetector equipment for the technique may be arranged to producenecessary measurement data for a reconstruction procedure utilizingappropriately-selected spatial parameters to thereby yield the desiredone millimeter spatial resolution. Finally, metabolic imaging with thetechnique should be realizable if imaging as localized spectroscopy isenvisioned in the sense that each point in the image is assigned anabsorption spectrum. Such an assignment may be used, for example, tomake a map of oxygenation by measuring the absorption spectra forhemoglobin at two different wavelengths, namely, a first wavelength atwhich hemoglobin is saturated, and a second wavelength at whichhemoglobin is desaturated. The difference of the measurements can yielda hemoglobin saturation map which can, in turn, give rise to tissueoxygenation information.

The first proposals for optical imaging suggested a mathematicalapproach (e.g., backprojection algorithm) that is similar to that usedto generate X-ray computerized tomography images. Light from a pulsedlaser is incident on the specimen at a source position and is detectedat a detector strategically placed at a point to receive transmittedphotons. It is assumed that the earliest arriving photons (the so-called"ballistic photons") travel in a straight line between the source anddetector, and the transmitted intensity is used in a mathematicalreconstruction algorithm. In effect, only unscattered incident waves areconsidered as being useful for forming an image of any objects embeddedin the specimen and, accordingly, techniques are employed to eliminatescattered light from the detection process, such as arranging a detectorwith "fast gating time" to only process the earliest arriving photons.However, since it is known that the ballistic photons are attenuatedexponentially, if the specimen has a thickness exceeding a predeterminedvalue, imaging is virtually impossible in many practical situations.

The latest proposals for optical imaging are now directed toward imagingsystems which use scattered and diffused radiation to reconstruct arepresentation of the interior of a specimen. Representative of priorart in this field is U.S. Pat. No. 5,070,455 issued to Singer et al(Singer) on Dec. 3, 1991. The system disclosed by Singer uses radiation,such as photons or other particles, which will be scattered to asignificant degree by the internal structure of a specimen. In thesystem, a specimen is irradiated and measurements of the attenuated andscattered radiation are effected at a number of points along theexterior of the specimen. It has been determined by Singer that suchmeasurements are sufficient to determine the scattering and attenuationproperties of the various regions inside the specimen. In accordancewith the disclosure of Singer, the interior of the specimen is modeledas an array of volume elements ("voxels"). Each voxel in the model ofthe specimen has scattering and attenuation properties which arerepresented by numerical parameters that can be mapped so as to generateseveral images of the interior of the specimen.

The particular technique used by Singer to reconstruct the interior ofthe specimen can best be characterized as an "iterative" procedure. Thisprocedure is now described in some detail so as to pinpoint itsshortcomings and deficiencies. After collecting the imaging data, thescattering and attenuation coefficients for the voxels are assignedinitial values, which helps to shorten the computation process--butwhich is also the characteristic of iterative or non-direct solution toa mathematical minimization problem. Next, the system computes theintensity of light that would emerge from the specimen if the interiorof the object were characterized by the currently assigned values forthe scattering and attenuation coefficients. Then, the differencebetween the measured light intensities and the computed lightintensities are used to compute an "error function" related to themagnitude of the errors of reconstruction. This error function (alsocalled "cost function" in minimization procedures) is then minimizedusing a multi-dimensional gradient descent methodology (such asFletcher-Powell minimization), i.e., the coefficients are modified so asto reduce the value of the error function.

The process of computing exiting light intensities based on thecurrently assigned values for the scattering and attenuationcoefficients, and then comparing the differences between the computedvalues and measured values to generate a new approximation of thescattering and attenuation properties of the interior of the specimen,continues until the error function falls below a specified threshold.The final values of the scattering and attenuation coefficients fromthis process are then mapped so as to generate a series of images of theinterior of the specimen, thereby depicting the attenuation andscattering characteristics of the specimen's interior--which presumablywill disclose both normal and abnormal conditions.

Singer thus discloses a technique to reconstruct an image by inversionusing an iterative minimization procedure. Such an approach is moreformally characterized as a "heuristic", in contrast to an "algorithm",since no verification or proof of even the existence of a solution usingthe approach has been offered. There are essentially an infinite numberof scattering and attenuation coefficients under such a regime, andthere is absolutely no assurance that the particular coefficientsdetermined using the iterative technique are the actual coefficients forthe specimen's interior. Moreover, such a heuristic method has a highcomputational complexity which is exponential in relation to the numberof voxels and which is, in turn, a characteristic of difficultoptimization problems with many local minima. The computationalcomplexity of such a approach renders the reconstruction methodvirtually useless for imaging.

There are other approaches presented in the prior art which are closelyrelated to that presented by Singer; these approaches also effect anindirect inversion of the forward scattering problem by an iterativetechnique which provide little, if any, physical insight.

Representative of another avenue of approach in the prior art is thesubject matter of U.S. Pat. No. 5,213,105 issued to Gratton et al(Gratton). With this approach, a continuous wave source of amplitudemodulated radiation irradiates an object under study, and radiationtransmitted or reflected by the object is detected at a plurality ofdetection locations, as by a television camera. The phase and theamplitude demodulation of the radiation is detected, and a relativephase image and a demodulation amplitude image of the object aregenerated from, respectively, the detected relative phase values anddetected demodulation amplitudes of the radiation at the plurality oflocations. However, while Gratton does generate data from a continuouswave source, there is no teaching or suggestion of a methodology to usethe data to directly reconstruct the image from the data--rather thedata of Gratton is merely used to obtain and display, separately andindependently, the shift in relative phase and the change in modulationof the wavefront of amplitude modulated radiation as a result ofpropagation of radiation through a scattering medium. Thus, there is noteaching or suggestion in Gratton of how the totality of the spatialmodulation and phase data can be combined in a unified approach todirectly reconstruct the image. It should be noted that at highmodulation frequencies the Gratton procedure is essentially equivalentto ballistic imaging, while at low modulation frequencies an entirelynew approach to the reconstruction problem is necessary.

SUMMARY OF THE INVENTION

These limitations and other shortcomings and deficiencies ofconventional techniques are obviated, in accordance with the presentinvention, by utilizing a direct reconstruction methodology, andconcomitant system, to generate an image of an object underinvestigation; the direct reconstruction formulation guarantees both theexistence and uniqueness of the imaging technique. Moreover, the directreconstruction method significantly reduces computational complexity.

In accordance with the broad aspect of the present invention, the objectunder study is irradiated by a continuous wave source at a givenfrequency and the transmitted intensity due predominantly to diffusivelyscattered radiation is measured at selected locations proximate to theobject wherein the transmitted intensity is related to the diffusioncoefficient by an integral operator. The image representative of theobject is directly reconstructed by executing a prescribed mathematicalalgorithm, determined with reference to the integral operator, on thetransmitted intensity measurements. In addition, radiation at differentwavelengths effects imaging as localized spectroscopy.

The organization and operation of this invention will be understood froma consideration of the detailed description of the illustrativeembodiment, which follows, when taken in conjunction with theaccompanying drawing.

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 exemplifies the transmission of light through a specimencontaining an absorbing object in the ballistic limit;

FIG. 2 exemplifies the transmission of light through a specimencontaining an absorbing object in the diffusion limit;

FIG. 3 depicts an object embedded in a medium for the case of withconstant absorption and fluctuating diffusion.

FIG. 4A-4C depict plots of the diffusion kernel for the characteristicintegral equation;

FIG. 5A and 5B illustrate an object in a medium and a simulatedreconstruction of the object using the direct reconstruction techniqueof the present invention;

FIG. 6 illustrates a high-level block diagram of an illustrativeembodiment of the tomography system in accordance with the presentinvention;

FIG. 7 is a high-level flow diagram of the methodology of the presentinvention;

FIG. 8 is a flow diagram depicting one methodology for computing thediffusion coefficient of the object under investigation; and

FIG. 9 shows a reconstructed exemplary object for diffusion imaging.

The same element appearing in more than one FIG. has the same referencenumeral.

DETAILED DESCRIPTION

To place in perspective the detailed description of the presentinvention and thereby highlight the departure from the art as disclosedand claimed herein, it is both instructive and informative to first gaina basic understanding of the imaging environment in which the presentinvention operates by presenting certain foundational principlespertaining to the subject matter in accordance with the presentinvention. Accordingly, the first part of the description focuses on ahigh-level discussion of the imaging systems relevant to the inventivesubject matter; this approach has the advantage of introducing notationand terminology which will aid in elucidating the various detailedaspects of the present invention. Abater this overview, the systemaspects of the present invention, as well as the concomitantmethodology, are presented with specificity.

Overview of the Present Invention

Multiple scattering of light presents a fundamental physical obstructionto optical imaging. The inventive subject matter of the presentinvention addresses this phenomena, with the surprising result thatdiffusive light contains sufficient information to image the opticaldiffusion coefficient of a highly scattering medium. This conclusionobtains from a version of inverse scattering theory that is applicableto multiple scattering in the diffusion limit. Using thisrepresentation, the first direct reconstruction procedure ever devisedfor imaging the optical diffusion coefficient of a highly scatteringmedium probed by diffusing waves is elucidated. In contrast totechniques which utilize unscattered (ballistic) photons for imageformation, the procedure in accordance with the present invention allowsfor the imaging of objects whose size is large compared to the averagescattering mean free path.

The familiar opaque or cloudy appearance of many objects havingimpinging light may be explained by the phenomenon of multiple lightscattering. (It is to be noted that terminology will be generalizedhereinafter so that an "object" is the physical manifestation of what isunder study--such an object may stand alone, may be embedded in aspecimen or a sample; in any case, the context of the descriptivematerial about an object will be set forth with clarity the meaning tobe attached to the generic term "object" in that context.) Thedisclosure and teachings of the present invention address the problem ofimaging an extended object that is embedded in a highly scatteringmedium. Since diffusively transmitted light contains sufficientinformation for direct image reconstruction, the problem can beexpressed in a tractable form amenable to an essentially closed-formsolution--meaning that there is no need to rely upon or resort to aniterative/minimization-type reconstruction with all its shortcomings andpitfalls.

To elucidate the direct reconstruction process at its most fundamentallevel, a simplified absorbing system to which direct reconstruction isapplicable is first described, namely, one in which a plane wave oflight (photons) of wavelength λ is incident upon a sample of lineardimension L that contains a spatially-extended absorbing objectcharacterized by a position-dependent optical absorption function; thewidth L is aligned with the impinging incident wave. If it is furtherassumed that photons are scattered by particles whose size is largecompared to λ, then the scattering is described by a transport mean freepath, designated l*; the mean free path characterizes the averagedistance a photon travels before its direction is randomized. In thesingle-scattering regime, that is, where l*>>L, it is observed that mostof the incident wave is unscattered upon exiting the sample and thus maybe used to form a projection image of the absorbing object; this effectis depicted in FIG. 1. In FIG. 1, light rays 101 of wavelength λ impingeon front 105 of sample 110 containing absorbing object 120, wherein thelight rays transmitted through sample 110 exiting back 106 of sample 110form a projection image represented by trace 130. The transmittedintensity represented by trace 130 is related to the line integral ofthe optical diffusion coefficient along the direction of propagation ofthe unscattered wave. This gives rise to the so-called Radon transformof the diffusion coefficient. By inversion of the Radon transform, it ispossible to recover the diffusion coefficient and thus an image ofabsorber 120 is reconstructed. As already alluded to above, allcommercially available imaging techniques are based on this simplephysical principle.

In the multiple-scattering regime, that is, where l*<<L, a wave scattersmany times while traversing the sample. In this situation, with λ<<*l,the path of a single photon may be described as a diffusive random walkwhere D=1/3(c/n)l* is a suitable diffusion constant, with c being thespeed of light, n being the index of refraction, and c/n being the speedof light in the medium of the sample. The unscattered, or ballisticphotons, are exponentially attenuated with a static transmissioncoefficient T_(ball) ˜exp(-L/l*). The dominant contribution to thetransmitted intensity is provided by diffusive photons with a diffusivetransmission coefficient T_(diff) ˜l*/L which, even with coherentillumination, forms a complicated interference pattern that does notcontain a simple image of the sample; such a pattern is illustrated inFIG. 2 (which has essentially the same pictorial representation as FIG.1, except that the physical system of FIG. 2 is such that l*<<L ascontrasted to l*>>L in FIG. 1). In FIG. 2, light rays 201 of wavelengthx impinge on front 205 of sample 210 and eventually exit sample 210 fromback 206. Absorbing object 220 gives rise to trace 230, which isrepresentative of the complicated transmitted light pattern exiting back206. In accordance with the present invention, there is devised aclosed-form procedure for utilizing the information in such complicatedpatterns as exhibited by trace 230 to locate an object and thus performoptical imaging in the multiple-scattering regime.

Indeed, it has frequently been pointed out in the prior art thatballistic photons convey the least distorted image information whilediffusive photons lose most of the image information. For this reasonseveral elegant experimental techniques have been designed to select theballistic photon contribution either by optical gating, holography, orfiltering of the diffusive photons by optical absorption. There is,however, an intrinsic physical limitation of any technique that reliessolely on ballistic photons. This may be appreciated by considering theexponential attenuation of ballistic photons relative to the mildalgebraic attenuation of diffusive photons. In particular, if the samplesize L is sufficiently large compared to l*, then T_(ball) will fallbelow an experimentally measurable threshold (e.g., if l* is about 0.5millimeters, then the attenuation is proportional to e⁻⁴⁰ in only 2centimeters).

Thus, the likelihood of now reconstructing important and valuable imagesheretofore believed to be virtually impossible to reconstruct provides astrong motivation to overcome the limitations of ballistic imaging byemploying multiply scattered diffusive photons for image reconstruction.From fundamental physical principles, such a reconstruction from theinterference pattern of diffusive transmitted light is attainable sincesuch reconstruction is uniquely determined by two parameters, namely,the absorption and diffusion coefficients of the highly scatteringsystem. As presented herein in accordance with the present invention,the diffusive transmission coefficient is related to the diffusioncoefficient by an integral equation. Then the image may be directlyreconstructed using a suitable algorithm which references this integraloperator. In contrast to ballistic methods, the resulting reconstructionmay be used to image samples whose size L is large compared to l*.

Function Theoretic Basis for Diffusion Imaging

Consider the propagation of a diffusing wave in a highly scatteringmedium due to an amplitude-modulated (AM), continuous wave (CW) source.As an example, such a source may be expressed as before, that is,S(x,t)=(1+Ae^(i)ωt)S(x), where S(x,t) is the source power density, ω isthe radian frequency of the continuous wave, and generally A<1.Typically, for ω=2πf, f is in the radio frequency (RF) range of 100 MHzto 1 GHz. For a point source at the origin, S(x)=S₀ δ(x) where S₀ is thesource power. Only the behavior of the energy density of the system isconsidered.

The energy density u(x,t) of a diffusing wave for constant absorptionand fluctuating diffusion obeys the following relation:

    ∂.sub.t u(x,t)=∇((D.sub.O +D(x))∇u(x,t))-α.sub.O u(x,t)+S(x,t),      (1)

where D₀ is the background diffusion of the medium, D(x) is thefluctuation in diffusion away from the background, and the absorptioncoefficient of the medium and the object is presumed to be constant (α₀)for this formulation. The diagram of FIG. 3 depicts these relations,namely, object 310 is shown immersed in medium 311 which has constantabsorption α₀ and diffusion constant D₀ ; object 310, on the other hand,has diffusion coefficient D(x) and constant absorption α₀. (Also shownfor completeness is i^(th) source 312 and j^(th) detector 313surrounding object 310, as discussed shortly).

Since only the long-time solution to the diffusion equation expressed inequation (1) is of concern, a solution of equation (1) in the followingform is expected:

    u(x,t)=u.sub.o (x)+Ae.sup.iωt u.sub.107 (x),         (2)

where u_(o) (x) is the zero frequency (DC) solution and u.sub.ω (x) isthe oscillatory (AC) solution at the angular frequency ω.

In a typical experiment, u_(ac) (x,t) is derived from measurementswhere:

    u.sub.ac (x,t)=e.sup.iωt u.sub.107 (x).              (3)

It is noted that u_(ac) (x,t)=|u.sub.ω (x)|e^(i) (¹⁰⁷ t +φ(x)) and that|u.sub.ω (x)| (the "modulus") and φ (x) (the "phase") are separatelymeasurable with appropriate detectors, as discussed shortly.

The relation expressed in equation (3) leads to the following definitionof transmission coefficient--defined as the transmitted intensitythrough the diffusing medium when the incident diffusing wave has unitflux:

    T.sub.ac =u.sub.ac (x)/u.sub.ac.sup.0 (x)=u.sub.ω (x)/u.sub.ω.sup.0 (x),                              (4)

where u.sub.ω⁰ (x) is the oscillatory part of the energy density whenthe diffusing object is not present. Thus, based on equation (4), it issufficient to determine u.sub.ω (x) to obtain T_(ac). Accordingly, itcan be shown in the weak fluctuation limit that ##EQU1## is designatedas the "diffusion kernel". Here, G₀ (x,x') is the unperturbed Green'sfunction for the diffusing wave, which in an infinite medium has thefollowing form: ##EQU2## Integral equation (5) relates the transmissioncoefficient of the diffusing wave to the diffusion coefficient.

An analytical expression for Γ_(D) (x;x₁,x₂) for an infinite medium maybe obtained with the result ##EQU3## Contour plots of Γ_(D) (x;x₁,x₂)are shown in FIG. 4 for various exemplary values of ωτ_(D), where τ_(D)=L² /D; in particular, FIGS. 4A-4C show Γ_(D) (x;x₁,x₂) for ωτ_(D) equalto 10⁴, 10², and 1, respectively. As depicted, at high frequencies (FIG.4A)) the diffusion kernel 401 is largely concentrated on the lineconnecting the source and detector. At lower frequencies (FIGS. 4B and4C), in the multiple scattering regime, the kernel (402 and 403,respectively) accounts for more extensive spatial contributions. Thusthe kernel provides a physical picture of photon transport in the CWapproach in the diffusion limit.

The central problem in photon diffusing wave imaging is thereconstruction of the image from transmission measurements for a familyof source-detector pairs. The description of a suitable reconstructionprocedure requires the solution of the fundamental integral equation(5). This integral equation is a Fredholm equation of the first kind.Such equations are typically ill-posed and it is well-known that theirsolution requires the introduction of a regularization method. Such aregularized solution of equation (5) may be obtained by singular valuedecomposition and is given by ##EQU4## is the regularized generalizedinverse of Γ_(D) (x;x₁, x₂), and where sources and detectors areexcluded from the volume of integration, that is, integration isperformed over a measurement surface surrounding the object. Here,σ_(n), f_(n), g_(n) denote the singular values and correspondingsingular functions of Γ_(D), Γ* _(D) *Γ_(D) f_(n) =σ_(n) ² f_(n), Γ_(D)f_(n) =σ_(ngn), β is a regularization parameter and R.sub.β is asuitable regularizer. Typically, R₆₂ is taken to be R.sub.β (σ)=σ/(β+σ²)so that small singular values are cut off smoothly. Equations (10) and(11) give the formal solution to the image reconstruction methodology indiffusion tomography.

The existence of the explicit inversion formula provided by equation(10) is of clear importance for the development of an imagereconstruction algorithm. The inversion formula, however, must beadapted so that transmission measurements from only a finite number ofsource-detector pairs may be used. One approach to this problem is toconsider a direct numerical implementation of the regularized singularvalue decomposition in equation (10). Here the integral equation (11) isconverted into a system of linear equations by an appropriatediscretization method such as collocation with piecewise constantfunctions. This method requires that measurements of the transmissioncoefficient be obtained from multiple source-detector pairs; each paircontributes multiple frequency points as well, such as by varying (oover the 100-500 MHz range in increments of 50 MHz. Thus at least asmany source-detector pair/frequency point combinations are required aspixels in the reconstructed image. It is important to appreciate thatthe computational complexity of such a real-space reconstructionalgorithm is O(N³) where N is the number of pixels in the reconstructedimage. It is noted that this is simply the complexity of the associatednumerical singular value decomposition.

The direct reconstruction of a two-dimensional object using computerprocessing on simulated transmission data is shown in FIG. 5; thereconstructed object of FIG. 5 serves as a pictorial visualization ofthe effectiveness of the direct reconstruction process. In FIG. 5A,object 520 is shown as embedded in medium 510. Since the shape of object520 is known, it is possible to mathematically describe, and thereby tocalculate the emanation of photons from back 506 due to photonsimpinging on front 505, that is, solve the so-called forward problem indiffusing wave imaging. Given the transmission intensity of photonsdetected proximate to back 506, the so-called inverse problem can besolved to directly reconstruct an image of object 520 --such a directlyreconstructed image is depicted by object 521 in FIG. 5B. It isimportant to emphasize that, although the transmission data issimulated, the algorithm used in the computer processing to reconstructimage 521 is the very one used to process actual measurements oftransmission intensity detected from an actual sample. Such simulationsafford the opportunity to study, for example, noise effects on thetransmission intensity data and sensitivity of the reconstructiontechnique to the various locations of the source-detector pair.

Details of the Diffusion Imaging Aspect of the Present Invention

A. SYSTEM

As depicted-in high-level block diagram form in FIG. 6, system 600 is adirect reconstruction imaging system for generating an image of anobject using measurements of transmitted radiation (e.g., photons)emanating from an object in response to photons impinging on the object.In particular, object 610 is shown as being under investigation. System600 is composed of: CW source 620 for irradiating object 610; dataacquisition detector 630 for generating the transmitted intensity ofradiation emanating from object 610 at one or more strategic locationsproximate to object 610, such transmitted intensity being determinedfrom measurements of both the modulus and phase (to obtain T_(ac)) via aphase-sensitive detector; position controller 640 for controlling thelocation of detector 630 relative to source 620; and computer processor650, having associated input device 660 (e.g. a keyboard) and outputdevice 670 (e.g., a graphical display terminal). Computer processor 650has as its inputs positional information from controller 640 and themeasured transmitted intensity from detector 630.

Source 620 and data acquisition detector 630 of system 600 are realizedwith conventional components; illustratively, Gratton discloses suchsource and detector arrangements, and the teachings of Gratton withrespect to the measurement procedure are hereby incorporated withreference. Thus, as before, CW source 620 is composed of a conventionalnear infrared laser operating in the 600 to 1200 nm region and at apower level of approximately 10 watts such as provided by a mode-lockedNeodymium YAG laser (Antares, model 76S-ML-SHG). Detector 630 is, forinstance, composed of a image intensifier (such as provided by PrincetonInstruments, Inc., OMA detector model IRY512 G/RB) which feeds aCCD-type television camera (such as provided by ITT, Fort Wayne, Ind.,model CCD F144). With the arrangements disclosed by Gratton, both thephase and modulus of the diffusively scattered radiation may begenerated.

Position controller 640 is utilized whenever CW source 620 and/or dataacquisition detector 630 may be composed of a plurality of radiationsources and detectors, respectively, in order to control which of theplurality of sources may be energized for a given time period and whichof the plurality of detectors may be activated during a prescribed timeinterval. As will be discussed in more detail below, in a practicalimplementation of the direct reconstruction imaging technique, it isoftentimes necessary to measure the transmitted intensity effected by anumber of source-detector positions surrounding object 610. For the sakeof expediency, generation of the required transmitted intensity data isexpeditiously accomplished by having arrays of P laser sources and Qphoton detectors. Accordingly, source 620 may be composed, in its mostgeneral implementation, of P CW sources or the like arrangedstrategically around the periphery of object 610. Similarly, dataacquisition detector may be composed, in its most general realization,of Q radiation detectors or the like also arranged strategically aroundthe periphery of object 610 and in a cooperative relation with the Psources.

The point of departure between the inventive subject matter herein andthe prior art resides in the processing of the measured data. Computer650 stores a computer program which implements the direct reconstructionalgorithm; in particular, the stored program processes the measuredtransmitted intensity data to produce the image of the object understudy using a prescribed mathematical algorithm; the algorithm isdetermined with reference to the integral operator relating thetransmitted intensity to the diffusion coefficient. The processingeffected by computer 650 is the focus of the discussion of themethodology section of this description, which follows immediately.

B. METHODOLOGY

B.1 Computational Model

The fundamental integral equation expressed by equation (5), repeatedhere,

    -InT.sub.ac (x.sub.1,x.sub.2)=∫d.sup.3 xΓ.sub.D (x;x.sub.1,x.sub.2)D(x),                                  (5)

is in the form of a Fredholm equation of the first kind (specificallyreferred to herein as the Schotland's Second Frequency-Domain IntegralEquation). Such equations are typically written in the form Kf=g, or

    ∫K(x,x')f(x')d.sup.3 x'=g(x)                          (12)

where f,g are elements of appropriately selected function spaces.Equation (12) is said to be ill-posed if (a) it is not solvable, (b) aunique solution does not exist, or (c) the solution does not dependcontinuously on the data. The latter case (c) is of primary interest inthe numerical study of ill-posed problems because it may lead tonumerical instability. This is particularly important if the data isimprecisely known or is the subject to statistical uncertainties, suchas measurement inaccuracy or noise, which would be the situation formeasurements for imaging. There are methods for conditioning ill-posedproblems. First, if the solution does not exist, the minimizer of||Kf-g|| is defined as a solution. Non-uniqueness is handled by choosingthe minimizer with the least norm. Finally, continuity is restored byintroducing "regularization" to the solution procedure.

Solving for the minimizer with the least norm yields the "normalequation" relating to equation (12); the normal equation is of the form

    K*Kf=K*g,                                                  (13)

where K* is the adjoint of K, and the property that K*K is self-adjointhas been employed. Thus, a solution for f in equation (12) is of thefollowing form:

    f=(K*K ).sup.-1 K* g ≡=K+g.                          (14)

From equation (14),

    K+=(K*K).sup.-1 K*                                         (15)

is called the "generalized inverse" of K.

B.2 Singular Value Decomposition

If K is such that a mapping from H₁ to H₂ occurs, where H₁ and H₂ areHilbert spaces, then K*K is a self-adjoint, positive operator. If theeigenfunctions and eigenvalues of K*K are denoted by {f_(n) } and {σ_(n)² }, respectively, then the following relation obtains: K*Kf_(n) =σ_(n)² f_(n).

The {σ_(n) } are the singular values of K. Also, the {f_(n) } form abasis for H₁. The singular values are ordered as σ₁ ² ≧σ₂ ² ≧ . . . ≧0,where multiplicities are counted and 0 can appear with infinitemultiplicity.

If {g_(n) } is defined by

    Kf.sub.n =σ.sub.ngn,                                 (16)

then the {g_(n) } are a basis for Hilbert space H₂. Moreover, it thenfollows that

    K* g.sub.n =σ.sub.n f.sub.n. (17)

To derive the singular value decomposition of K, put K in the form

    K=I.sub.H2 KI.sub.H1                                       (18)

and use the identities, ##EQU5## where x denotes the tensor product.Manipulation of equations (18)-(20) leads to ##EQU6## Equation (21) iscalled the "singular value decomposition" of K.

The singular value decomposition of equation (21) can now be used toobtain a form for the generalized inverse K+ of equation (15). As aresult of equation (21), ##EQU7## then it directly follows, aftersubstitution of equations (22) and (23) into equation (15), that##EQU8## Now, using equations (14) and (24), the solution of Kf=g isf=K+g, which is of the form ##EQU9##

If some of the σ_(n) 's vanish, then K+ is not well-defined and, inparticular, is not continuous. To resolve this anomaly, theregularization procedure is introduced.

B.3 Regularization

To condition the singular value decomposition, the following expressionis now defined: ##EQU10## where the regularizer R.sub.β (σ) has theproperties

    (i) R.sub.62 (σ)=1/σ as β→0.sup.+ ;

    (ii) R.sub.β (σ)˜1/σ for σ>>0 (with β>0);(27)

    (iii) R.sub.β (σ)→0 as σ→0(with β>0).

For instance, two natural choices (others are possible) include:

    (a) R.sub.β (σ)=1/σ for σ>β; otherwise, R.sub.β (σ)=0;                                 (28)

    (b) R.sub.β (σ)=σ/(β+σ.sup.2). (29)

(One typical heuristic criterion is to set β˜O(σ₁)).

Thus the solution of equation (12) may be written as ##EQU11## where

    K.sub.β.sup.+ (x,x')=Σ.sub.n R.sub.β (σ.sub.n)f.sub.n (x)g.sub.n (x').                                          (31)

(The form of equation (11) follows from the generic notation used toobtain equation (31)).

B.4 Numerical Solution of the Schotland's Second Integral Equation

The above developments for the formal solution of a general Fredholmequation of the first kind, including the techniques of singular valuedecomposition and regularization, may now be applied to implement thenumerical solution of equation (11)--Schotland's Second IntegralEquation:

    -InT.sub.ac (x.sub.1,x.sub.2)=∫d.sup.3 xΓ.sub.D (x;x.sub.1,x.sub.2)D(x),                                  (5)

the formal solution, by way of summary, is given by equations (10) and(11), as follows: ##EQU12##

For a three-dimensional object, denoted Ω, it is supposed that there areP sources and Q detectors used to probe the object. These sources arespaced about the periphery of the object and, operating in conjunctionwith the sources, there are suitably placed detectors. For the sake ofsimplicity, a single frequency is considered in the followingexposition. In general, the results may be readily extended to the caseof multiple frequencies. Let i, i=1,2, . . . ,P and j, j=1,2, . . . ,Qbe indices corresponding to the P sources and Q detectors; then, for agiven frequency, equation (5) becomes: ##EQU13## (In equation (32) andfor the remainder of this section, T≡T_(ac) and Γ≡Γ_(D).) Now Γ and Dare discretized by decomposing Ω into "voxels" (i.e., volume elementshaving basically equal sides) B_(m), m=1,2, . . . , M which cover theobject. It is then assumed that the granularity is such that D and Γ areconstant in each box. To recast equation (32) in a standard form, thefollowing identifications are made:

    |B.sub.m |Γ.sub.ij.sup.t (x.sub.m)≡A.sub.ij.sup.m,                           (33)

    D(x.sub.m)≡a.sub.m,                                  (34)

and

    -InT.sub.ij.sup.t ≡b.sub.ij,                         (35)

where |B_(m) |is the volume of a voxel, and a_(m) is the strength ofD(x_(m)) at the middle of the m^(th) voxel. Then, using thesedefinitions, equation (32) becomes ##EQU14## for m=1,2,. . . ,M; i=1,2,. . . ,P; and j=1,2,. . . ,Q. In matrix form, equation (36) isrepresented as Aa=b, where A is a (PQ by M) matrix, so equation (36)gives PQ equations in M unknowns. Thus, there must be at least as manysource-detector pairs (PQ) as voxels (M). It is preferable to"overdetermine" equation (36) by having PQ>M, or by using multiplefrequencies for each source-detector pair. If there are K frequencies,then matrix A is a (KPQ by M) matrix. Typically in practice KPQ=3M.There are many possible ways to arrange for desired number of frequencymeasurements. For instance, it is possible to fix ω, and arrange formultiple sources and detectors each operating at ω, but only operate onesource at a time. As another example, it is possible to vary ω, andarrange for multiple sources and detectors, each detector being arrangedto detect the changes in ω(e.g., by selecting one of a plurality ofband-pass filters), but only operate one source at a time. Finally, itis possible to employ multiple frequencies emitted by multiple sources,and a single detector tuned simultaneously to each of the multiplefrequencies, with all sources operating simultaneously.

The solution of singular value decomposition applied to a matrixformulation is a well-known technique. For example, as previouslyindicated, a procedure for singular value decomposition is described inthe text "Numerical Recipes", by Press, Flannery, Teukolsky, andVettering, 1986, Cambridge University Press, Cambridge, England. Acommercially available software package implementing the singular valuedecomposition, called Interactive Data Language (IDL) available fromResearch Systems Inc. of Denver, Colo., may be used in practice; IDL wasspecifically designed for scientific computations, especially imageprocessing applications. With IDL, a subroutine-like call of the form"SVD Matrix!" (e.g., SVD A! in terms of the above A matrix) returns thesingular values as well as the quantities, denoted the projectionoperators, from which {f_(n) } and {g_(n) } obtain.

Once the singular value decomposition has been effected, regularizationaccording to equation (33) is readily accomplished in order to obtainthe regularized, generalized inverse which, for the matrix A, is denotedA+. The solution to the discretized Schotland's Second Integral Equationbecomes a=A+b.

c. FLOW DIAGRAMS

The methodology discussed in the previous section is set forth inhigh-level flow diagram 700 in FIG. 7 in terms of the illustrativeembodiment of the system shown in FIG. 6. With reference to FIG. 7, theprocessing effected by control block 710 enables photon source 620 anddata acquisition system 630 so as to measure energy emanating fromobject 610 due to CW source 620. These measurements are passed tocomputer processor 650 from acquisition system 630 via bus 631. Next,processing block 720 is invoked to retrieve the pre-computed and storeddiffusion kernel as expressed by equation (6). In turn, processing block730 is operated to execute the direct reconstruction algorithm set forthwith respect to equations (12)-(36), thereby determining the diffusioncoefficient D(x). Finally, as depicted by processing block 740, thereconstructed image corresponding to the diffusion coefficient D(x) isprovided to output device 670 in a form determined by the user; device670 may be, for example, a display monitor or a more sophisticatedthree-dimensional video display device.

One illustrative manner of carrying out the direct reconstructionexhibited by block 720 is further depicted by high-level flow diagram800 of FIG. 8. In particular, processing block 810 shows that the firststep is to form the kernel matrix (the A matrix as determined bydiscretization, that is, A_(ij) ^(m) for m=1,2, . . . , M; i=1,2, . . .,P; j=1,2, . . . ,Q, of equation (33)). Next, processing block 820 isinvoked to compute the singular value decomposition of the kernel matrixA. Then, processing block 830 is executed to generate the regularized,generalized inverse A+. Finally, block 840 is invoked to obtain thesolution a=A+b, where a represents the discretized values of thediffusion coefficient.

With reference to FIG. 9, there is shown the reconstruction of anexemplary object. In particular, FIG. 9(a) shows an originaltwo-dimensional object 902 embedded in a 10 cm×10 cm specimen 901. FIGS.9(b)-(d) show the direct reconstruction of the image 902 (e.g., ,9021-9023 in FIGS. 9(b)-9(d), respectively) in the presence of additiveGaussian noise of 0.1%, 1%, and 5% to indicate the relativeinsensitivity of the direct reconstruction to noise. In this example, α₀=1.0 cm⁻¹ and the contour levels 903, 904, and 905 refer to thefluctuation in the optical diffusion measured in units of cm² ns⁻¹. Thetransmission coefficients were obtained using Monte Carlo simulations.

The system and methodology described utilizes the free-space model ofthe diffusion kernel (equation (9)) so that the kernel is pre-computedand stored in computer processor 650 for recall during thereconstruction process. This is appropriate when object 610 issurrounded by an arrangement, such as a thin, rubber-like containerfilled with a substance (e.g., , the commercially available medicalproduct called Intralipid), so that the arrangement provides a spatialextent external to the object that effectively gives rise to afree-space condition surrounding the object. The object's actualboundary (e.g., a human skull during imaging of brain) becomes merelyanother shape that is determined by the direct reconstruction procedure.Intralipid is useful because it is a colloidal substance whereinparticles in the range of 0.5 microns to 2 microns are suspended, andthe substance, as packaged, does not rapidly deteriorate; moreover, thel* of such a substance is readily measurable.

Whereas the above formulation has focused on the case of amplitudemodulation, that is, to ω≠0, the formulation also applies for so-calledDC imaging, that is, the case wherein to ω=0. Thus both cases, namely,AC imaging (ω≠0) and DC imaging (ω=0) are encompassed by the terminology"continuous wave."

Also, whereas the above formulation has focused on the case ofreconstruction from measurements of T_(ac), it is possible to formulatethe direct reconstruction methodology in terms of the transmissioncoefficient T_(mod) where

    T.sub.mod =|u.sub.ω (x)|/|u.sub.ω.sup.0 (x)|.(37)

With this formulation, only measurements of the modulus are required,thereby precluding reconstruction perturbations that may occur becauseof phase measurement inaccuracies in generating T_(ac). The integralequation using the modulus alone may be obtained as follows: ##EQU15##

Similarly, reconstruction based only on phase measurements may beperformed; in this situation, the following integral equation obtains:##EQU16## and where φ₀ (x₁,x₂) is the reference phase measured in theabsence of absorption.

Whereas the above discussion has focused on obtaining modulus and phaseinformation utilizing frequency domain measurements, it will be readilyappreciated by one of ordinary skill in the art that a time domainsource may be utilized to obtain, for example, a transmitted intensitydue to a pulsed source, and that such measured intensity may then beconverted to the frequency domain by a transformation process such as aFourier Transform.

It is to be understood that the above-described embodiment is simplyillustrative of the application of the principles in accordance with thepresent invention. Other embodiments may be readily devised by thoseskilled in the art which may embody the principles in spirit and scope.Thus, it is to be further understood that the methodology describedherein is not limited to the specific forms shown by way ofillustration, but may assume other embodiments limited only by the scopeof the appended claims.

What is claimed is:
 1. A method for generating a diffusion image of anobject having a variable diffusion coefficient, the method comprisingthe steps of:irradiating the object with a continuous wave source ofradiation, measuring a transmitted intensity due predominantly todiffusively scattered radiation wherein said transmitted intensity isrelated to the diffusion coefficient by an integral operator, anddirectly reconstructing the image by executing a prescribed mathematicalalgorithm, determined with reference to said integral operator, on saidtransmitted intensity; wherein said step of directly reconstructing theimage includes the step of computing a diffusion kernel.
 2. The methodas recited in claim 1 wherein the step of irradiating the objectincludes the step of successively irradiating the object with differentwavelengths.
 3. A system for generating a diffusion image of an objecthaving a variable diffusion coefficient, the systemcomprising:continuous wave radiation source means for irradiating theobject, detector means for measuring a transmitted intensity duepredominantly to diffusively scattered radiation wherein saidtransmitted intensity is related to the diffusion coefficient by anintegral operator, and means for directly reconstructing the image byexecuting a prescribed mathematical algorithm, determined with referenceto said integral operator, on said transmitted intensity; wherein saidmeans for directly reconstructing the image includes means for computinga diffusion kernel.
 4. The system as recited in claim 3 wherein saidsource means for irradiating the object includes means for successivelyirradiating the object with different wavelengths.